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Objectives  Our  overall  goal  is  to  develop,  implement  and  transfer  accurate  new 
numerical  methods  for  solving  moving  boundary  problems  in  materials  science.  We 
have  three  specific  objectives: 

•  Combine  semi-Lagrangian  time  stepping,  accurate  contouring  and  fast  geomet¬ 
ric  algorithms  to  develop  and  implement  accurate,  efficient  and  general  new 
methods  for  moving  sharp  interfaces. 

•  Develop  a  fast  modular  moving  interface  code  for  transfer  to  other  researchers, 
labs,  and  industry. 

•  Build  efficient,  accurate  and  general  integral  solvers  for  coupled  nonlinear  sys¬ 
tems  of  Poisson  and  heat  equations  modeling  common  material  phenomena,  and 
couple  these  solvers  to  our  modular  moving  interface  code. 

Status  of  Effort  We  have  attained  our  first  objective — an  accurate,  efficient  and 
general  moving  interface  method — with  the  semi-Lagrangian  contouring  method  re¬ 
ported  in  Publication  [PI],  We  are  developing  a  portable  C  code  with  fast  accurate 
contouring  techniques  based  on  deferred  correction  schemes,  and  fast  PDE  solvers 
based  on  integral  equation  techniques,  to  complete  our  second  and  third  objectives. 

Accomplishments  Our  work  on  moving  interface  problems  in  materials  science 
combines  fast  PDE  solvers  such  as  boundary  integral  methods  with  fast  geometric 
algorithms  and  semi-Lagrangian  implicit  representations  to  build  effective  new  nu¬ 
merical  methods. 

We  developed  an  implicit  boundary  integral  method  for  computing  periodic  den¬ 
drite  formation  in  the  symmetric  model  of  unstable  solidification  [7]  and  fast  algo¬ 
rithms  for  evaluating  heat  potentials  [2]  which  speeded  up  our  method  by  several 
orders  of  magnitude.  In  [6],  we  combined  the  boundary  integral  method  of  [7]  with 
fast  algorithms  from  [2,3,8]  and  the  level  set  method  of  [4]:  the  level  set  method  han¬ 
dled  topological  changes  effectively  while  fast  boundary  integral  techniques  ensured 


accuracy  and  efficiency  in  the  velocity  evaluation.  We  developed  and  analyzed  effi¬ 
cient  and  accurate  new  vortex  methods  for  modeling  convection  in  the  melt  [4,10,11], 
together  with  new  error  analyses  [12]  and  quadrature  rules  [9]  for  general  integral 
equations.  Since  1999,  we  have  focused  on  the  development  and  implementation  of 
highly  effective  new  numerical  methods  for  general  moving  interface  problems  and 
widely  applicable  subsidiary  computations.  We  summarize  three  projects  below:  the 
fast  modular  semi-Lagrangian  method  for  general  moving  interfaces  (described  in 
Publications  [P1,P2]  and  references  [13-16]),  accurate  contouring  methods  in  two  and 
three  dimensions,  and  fast  solution  of  two-point  boundary  value  problems  [P3,P4]. 


Moving  interface  problems  A  moving  interface  is  a  collection  T(t)  of  noninter¬ 
secting  oriented  closed  curves  (in  R2)  or  surfaces  (in  R3).  A  sufficiently  smooth 
moving  interface  has  an  outward  unit  normal  vector  N,  a  mean  curvature  C ,  and  a 
normal  velocity  vector  VN.  A  moving  interface  problem  specifies  VN  as  a  functional 
of  r(i).  Examples  include  passive  transport  V  =  N  ■  F  where  F(x,t )  is  given,  geo¬ 
metric  motion  V  =  p{&)  —  q{0)C  where  cos  6  =  N  •  x,  Ostwald  ripening  V  =  [J77  , 
where  u  solves  the  Laplace  equation  A u  =  0  off  r(f)  and  u  =  C  on  T(t),  and  models 
for  crystal  growth  V  =  [|^] ,  where  u  solves  the  heat  equation  ut  =  A u  off  T(t)  with 
the  geometric  boundary  condition  u  =  f(N ,  C,  V)  on  T(t). 


The  implicit  approach  Any  moving  interface  problem  can  be  reformulated  as  a 

PDE  for  a  function  <p  whose  zero  set  is  F(t).  The  normal,  curvature  and  velocity  are 

then  _  ~ 

Vp  „  „  Ar  t /  a r  <PtV<P 


N  = 


C  =  —V  •  N,  VN  = 


livin’  ~  liv^ll2' 

Given  an  extension  of  FAT  off  T(t)  to  a  globally  defined  function  W,  we  can  regard 
the  VN  formula  as  a  PDE  for  tp: 


<Pt  —  W  ■  =  0. 

We  solve  this  equation  on  an  adaptive  quadtree  mesh  to  eliminate  the  cost  of  going 
up  a  dimension.  Correct  viscosity  solutions  are  obtained  by  semi-Lagrangian  time 
stepping  with  exact  distancing  and  large  time  steps.  A  general  problem-independent 
velocity  extension  makes  our  method  modular  and  easy  to  apply. 

Semi-Lagrangian  methods  The  semi-Lagrangian  “CIR”  method  [1]  solves  cpt  — 
F(x,  t )  •  Vy  =  0  by  the  following  algorithm:  at  each  x  in  the  grid,  move  x  back 
with  velocity  F(x,tn )  to  s  =  x  +  kF(x,tn)‘,  interpolate  <p(x,tn)  to  the  point  s;  set 
<p(x,tn+ 1)  equal  to  the  interpolated  value.  Our  second-order  time  stepping  scheme 
couples  a  CIR  predictor  with  a  trapezoidal  corrector  using  the  velocity  evaluated 
from  the  CIR  approximation.  It  combines  the  unconditional  stability  of  CIR  with 
the  dramatically  reduced  dissipation  of  the  trapezoidal  rule.  Interpolation  error  is 
eliminated  by  exact  distance  finding  in  a  dynamic  quadtree  data  structure. 

We  first  tested  the  semi-Lagrangian  approach  in  [13]  by  solving  level  set  equations 
with  a  fixed  uniform  mesh  and  ENO  differencing.  Second,  we  developed  an  efficient 


new  redistancing  technique  with  the  aid  of  a  new  quadtree  structure  whose  cells  know 
the  signed  distance  to  nearby  elements  of  T(t).  The  quadtree  is  built  in  0(N  log  N) 
work  by  a  three-step  recursive  search  procedure.  Computations  of  moderate  complex¬ 
ity  are  speeded  up  400  times,  while  redistancing  the  CIR  method  on  a  uniform  mesh 
costs  considerably  less  than  moving  <p  one  step  [14].  Third,  the  semi-Lagrangian 
method  of  [15]  combined  CIR  with  an  adaptive  quadtree  mesh  to  build  a  method 
which  is  not  only  accurate  and  robust  but  also  optimally  efficient:  an  interface  T(£) 
with  N  degrees  of  freedom  costs  only  0(N log  N)  to  move  one  step.  Efficiency  is 
further  enhanced  by  the  unconditional  stability  of  the  semi-Lagrangian  time  stepping 
scheme:  large  time  steps  k  =  0(h)  can  be  taken  even  on  refined  meshes.  The  tree 
mesh  is  refined  with  a  functional  approach:  given  a  signed  distance  function  ip(x,tn), 
we  build  a  tree  at  time  £n+i  =  tn  +  k  by  recursive  evaluation  of  <p(x,  tn+ 1)  =  ip(s,  tn) 
at  projected  points  s  =  x  +  kF(x,tn).  The  criterion  for  splitting  a  tree  cell  is  simple: 
the  values  of  ip(x,  tn+ 1)  on  the  cell  are  smaller  than  the  size  of  the  cell. 

Our  semi-Lagrangian  method  for  moving  interfaces  [16]  combines  efficient  exact 
quadtree-based  redistancing,  stable  second-order  semi-Lagrangian  time  stepping,  a 
modular  problem-independent  velocity  extension,  and  exact  </?  interpolation  in  the 
CIR  scheme.  The  velocity  extension  technique  evaluates  the  nearest-point  extension 
on  a  distance  tree,  builds  a  continuous  interpolant,  and  satisfies  a  maximum  principle. 
Our  method  resolves  and  moves  complex  interfaces  at  optimal  cost  with  time  steps 
unconstrained  by  numerical  stability.  It  is  a  “black-box”  method  for  moving  inter¬ 
faces,  which  accepts  the  interface  and  its  velocity  at  time  t  and  returns  the  evolved 
interface  one  time  step  later.  Such  methods  simplify  moving  interfaces,  because  the 
numerics  are  independent  of  the  physical  problem  driving  the  interfacial  motion.  Nu¬ 
merical  results  show  that  the  method  converges  to  correct  viscosity  solutions  even  for 
difficult  moving  interface  problems  involving  merging,  faceting,  transport,  nonlocality 
and  anisotropic  curvature-dependent  geometry. 

A  fast  semi-Lagrangian  contouring  method  [PI]  General  moving  interface 
problems  are  solved  in  [PI]  by  a  new  approach:  extract  the  moving  interface  from 
an  explicit  semi-Lagrangian  advection  formula  with  efficient  geometric  algorithms 
and  fast  accurate  contouring  techniques.  A  modular  adaptive  implementation  with 
fast  new  geometry  modules  computes  highly  accurate  solutions  to  moving  interface 
problems  involving  merging,  anisotropy,  faceting,  curvature,  dynamic  topology  and 
nonlocal  interactions  of  PDE  type.  Exact  geometric  algorithms  are  tuned  for  speed; 
velocity  evaluation  and  time  stepping  are  efficiently  decoupled  from  interface  reso¬ 
lution;  fast  new  contouring  techniques  dramatically  increase  overall  accuracy.  An 
efficient  adaptive  framework  combines  the  high  resolution  of  front  tracking  with  the 
topological  robustness  of  implicit  representations.  Currently  three  projects  are  un¬ 
derway:  application  to  nonlocal  problems  such  as  Ostwald  ripening,  development  of 
accurate  general  contouring  schemes  in  two  and  three  dimensions,  and  fast  high-order 
solution  of  two-point  boundary  value  problems. 


Ostwald  ripening  [P2]  We  are  building  fast  nonlocal  velocity  evaluation  modules 
for  the  standard  moving  interface  problems  of  materials  science.  The  simplest  example 
is  Ostwald  ripening,  which  models  the  growth  of  larger  solid  drops  by  evaporation 
from  smaller  drops  with  total  solid  volume  conserved.  The  velocity  V  is  the  normal 
derivative  of  the  function  u  which  is  harmonic  off  T(t)  and  equal  to  the  curvature 

on  r (t).  We  evaluate  this  nonlocal  velocity  by  solving  the  integral  equation  of  classical 
potential  theory  and  applying  the  Dirichlet  to  Neumann  operator.  The  solution  u  is 
a  double  layer  potential 

of  an  unknown  density  (jl  on  T  =  T (i),  with  K(x,y)  =  ^log||a;  —  y||  the  free-space 
Green  function  for  the  two-dimensional  Laplace  equation.  The  density  y  solves  the 
integral  equation 

lM+Ld-mrMdv=°(x)’  x€r- 

Once  y  is  found,  it  is  convenient  to  view  the  harmonic  function  u  as  the  real  part  of  an 
analytic  function  U.  The  Cauchy-Riemann  equations  then  yield  V  as  the  tangential 
derivative  of  the  imaginary  part  of  U,  which  is  easier  to  compute  than  the  normal 
derivative  of  the  real  part  u.  Discretization  of  this  formulation  is  highly  accurate  if 
the  interface  is  represented  by  equidistant  points  in  arclength.  Detailed  resolution 
of  the  interface  requires  many  points,  so  fast  algorithms  such  as  the  fast  multipole 
method  play  an  important  role. 

Accurate  contouring  The  general  problem  of  finding  a  smooth  geometrically  con¬ 
strained  approximate  zero  set  of  a  function  which  can  be  evaluated  at  arbitrary  points 
occurs  frequently  in  computational  science  and  requires  a  robust  general  contouring 
package.  An  ideal  contouring  package  would  accept  function  values  (and  derivatives 
if  available)  at  arbitrary  points  and  produce  a  piecewise-smooth  approximation  to 
the  zero  set  with  corners  where  necessary.  Geometric  constraints  such  as  bounds 
on  curvature  away  from  corners  are  vital  in  applications  such  as  computer-controlled 
machining,  and  pose  a  major  complication  for  existing  public-domain  contouring  soft¬ 
ware.  We  are  developing  a  generally  useful  package  for  constrained  piecewise-smooth 
contouring  of  scattered  data  in  two  and  three  dimensions,  based  on  new  methods  of 
scattered  data  interpolation  and  a  new  floating-point  stability  analysis  of  the  con¬ 
touring  problem. 

Deferred  correction  solvers  for  boundary  value  problems  [P3,P4]  Contour¬ 
ing  on  a  subdomain  requires  the  solution  of  a  two-point  boundary  value  problem  for 
the  curve  connecting  known  zeroes  on  the  subdomain  boundary.  Thus  we  have  de¬ 
veloped  efficient  and  accurate  new  deferred  correction  techniques  for  the  solution  of 
the  two-point  boundary  value  problem 

y'  =  C(t)y  +  f(t),  a<t<b 

Ay(a)  +  By(b )  =  g 


for  a  vector-valued  function  y  :  [a,  b]  — »  Rq.  Deferred  correction  is  a  strategy  of 
systematically  promoting  a  low-order  scheme  such  as  the  midpoint  rule  to  an  efficient 
high-order  scheme,  by  applying  the  same  basic  scheme  to  solve  an  equation  for  the 
error  in  terms  of  the  residual.  The  error  equation  is  solved  efficiently  by  repeated 
use  of  a  highly  stable  block  arrowhead  QR  factorization.  This  approach  yields  stable, 
efficient  and  highly  accurate  schemes  of  orders  up  to  20,  with  naturally  adapted  grids. 
Fig.  1  illustrates  the  speed  and  accuracy  of  these  schemes,  applied  to  difficult  and 
singular  problems  which  challenge  standard  packages  [P3].  A  new  convergence  theory 
confirms  these  experimental  results  [P4]. 

The  future  Planned  future  developments  include  higher-order  accurate  contouring 
and  boundary  representations  with  constrained  geometry,  three-dimensional  imple¬ 
mentation,  and  applications  to  PDE-type  problems  where  the  moving  interface  is 
coupled  to  complex  materials  science. 
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